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ABSTRACT 

We investigated the orbital evolution of satellite galaxies using numerical sim- 
ulations. It has been long believed that the orbit suffers circularization due to the 
dynamical friction from the galactic halo during orbital decay. This circulariza- 
tion was confirmed by numerous simulations where dynamical friction is added 
as external force. However, some of the resent iV-body simulations demonstrated 
that circularization is much slower than expected from approximate calculations. 
We found that the dominant reason for this discrepancy is the assumption that 
Coulomb logarithm log A is constant, which has been used in practically all recent 
calculations. Since the size of the satellite is relatively large, accurate determina- 
tion of the outer cutoff radius is crucial to obtain good estimate for the dynamical 
friction. An excellent agreement between A^-body simulations and approximate 
calculations was observed when the outer cutoff radius is taken to be the distance 
of the satellite to the center of the galaxy. When satellite is at the perigalacticon, 
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the distance to the center is smaller and therefore log A becomes smaller. As a 
result, the dynamical friction becomes less effective. 

We apply our result to the Large Magellanic Cloud. We found that the 
expected lifetime of the LMC is twice as long as that would be predicted with 
previous calculations. Previous study predicts that the LMC will merge into the 
Milky Way after 7 G years, while we found that the merging will take place after 
14 G years from now. Our result suggests that generally satellites formed around 
a galaxy have longer lifetime than previous estimates. 

Subject headings: celestial mechanics, stellar dynamics — Galaxy:kinematics 
and dynamics — galaxies: Magellanic Clouds — Local Group — methods: nu- 
merical 



1. Introduction 

Recent observations have revealed that there are many satellite galaxies around the 
Milky Way. In the hierarchical clustering scenario, it is expected many of such dwarf satel- 
lites are formed. In fact, one of the most serious problems with the present hierarchical 
clustering scenario is that it predict too many satellite galaxies, about a factor of 10 more 
than the number observed in the Local group {e.g., Moore et ai, 1999). A number of expla- 
nations, including exotic theories which relies on hot or self-interacting dark matter, have 
been proposed. 

In this paper, we go back to the basic problem: how long are the satellites lives ? In other 
words, how do the orbits of satellites evolve through interaction with the gravitational field of 
its parent galaxy? The dominant driving force of the evolution is the dynamical friction. For 
satellites like the LMC-SMC pair and the Sagittarius dwarf, there are many detailed studies 
of their orbital evolution, in which the dynamical friction is included as the external force 
operating on the center-of-mass motion of the satellite. Well known works include Murai and 
Fujimoto (MCs, 1980) and Ibata and Lewis (Sagittarius, 1998). In both of these studies, and 
in all other studies where the dynamical friction formula is used, significant circularization 
of the orbit of the satellite was observed. This circularization is the natural result of the fact 
that the dynamical friction is proportional to the local density of the background stars, and 
therefore the strongest at the perigalacticon. 

However, recent A^-body simulations of the orbital evolution of satellites resulted in 
rather counter-intuitive result. Van den Bosch et al (1999, hereafter BLLS) performed the 
A^-body simulation of the satellite, where the parent galaxy was modeled directly as self- 
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consistent iV-body system. The satellite is modeled as one massive particle with spline 
potential softening used in PKDGRAV (Dikaiakos & Stadel, 1996). They investigated the 
evolution of the orbit for wide variety of model parameters such as the mass of the satellite 
and initial orbital eccentricity. They observed practically no circularization in any of their 
simulations. 

Jiang and Binney (2000, hereafter JB) performed fully self-consistent simulation of the 
satellite, where both the parent galaxy and the satellite are modeled as self-consistent TV- 
body systems. They compared their result with the result of approximate model in which 
the usual dynamical friction formula is used. Though they argued that the agreement is 
good, from their figure 3 it is clear that approximate models suffer stronger circularization 
and evolve faster than their iV-body counterpart. 

Neither of above two papers discussed the reason of this rather serious discrepancy 
between the result of iV-body simulations and previous analytic prediction. The purpose 
of this paper is to understand its cause. In section 2, we describe our model experiment 
designed to reproduce the discrepancy observed by BLLS and JB. In section 3 we show our 
result. Our result is consistent with both of the previous works. iV-body simulation showed 
only marginal circularization but approximate calculation using dynamical friction formula 
showed strong circularization. In section 4, we investigate the reason. There are several 
possible candidates for the reason. We consider a few of them, and found that a simple 
modification of the conventional form of the dynamical friction formula results in a quite 
remarkable improvement of the agreement between iV-body and approximate calculations. 
In section 5 we apply our formalism to the LMC. We found that the orbital evolution becomes 
significantly slower than prediction by previous calculations using conventional formula. For 
example, the lifetime of LMC was 7 Gyr with conventional formula, but is 14 Gyr with our 
formalism. We also discuss the implication of our result to the so-called "dwarf problem" . 

2. Numerical Simulation 

We carried out a set of numerical simulations to see whether the results obtained by 
BLLS and JB are really true or not. In this section, we describe the models we used. 

2.1. iV-body simulation 

We performed iV-body simulations of the evolution of a satellite orbiting in a massive 
dark halo of a galaxy. 
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The massive halo is composed of N equal mass particles, while the satellite dwarf is 
modeled by a single particle with a certain softening length. The softening is used to mimic 
the finite size of the satellite. 

We adopted a King model of the concentration ratio $ = 9 as a model of the galac- 
tic halo. The system of units is the Heggie unit (Heggie and Mathieu (1986)) where the 
gravitational constant G is 1, the mass and the binding energy are 1 and 0.25, respectively. 

JB used a composite disk+halo model in which the halo is expressed by particles and 
the disk is assumed to be rigid. BLLS used a single spherical halo. In both works, the halo 
density profile has the form 

r 2 c exp[-{r/r t ) k ] 

P = P° ^2^2 > C 1 ) 

where r c and r t are the core radius and the outer scale radius of the halo and po is the central 
density of the halo. BLLS adopted k — 2 while JB adopted k — 1. 

We did not follow the models. The standard dynamical friction formula is derived for 
the case of field stars with the Maxwell distribution. However, the distribution function 
associated with eq. (1) is rather different from the Maxwell distribution. This may cause 
difference in the effect of the dynamical friction. Also, the distribution function would 
relax to the Maxwellian through two-body relaxation, causing a small change in both the 
distribution function and the density profile. 

In addition, the range of radius for which the density slope is approximately —2 is rather 
narrow with this model, since the slope is noticeably shallower than —2 for r < 10r c . 

The distribution function of the King model is a simple lowered Maxwellian. Therefore 
the agreement with the true Maxwellian is very good within the half-mass radius. Also, 
since the distribution function is practically as close as the true Maxwellian as we can make, 
thermal relaxation is minimized, though it still presents (see e.g., Quinlan 1996). Also, the 
King model with \l>o = 9 has fairly wide range of radius in which the slope of the density is 
approximately —2. So it is a fairly good model for a spherical halo with flat rotation. 

The satellite galaxy is modeled by a single particle with mass M s and softening length 
e s . The force on the satellite from a particle in the halo is calculated as follows 

F CmM s (r sat - r halo ) 

(|r sat -r halo p + e 2 + e 2 aio) 3/2- U 

Here thaio is the softening length for the particle in the galactic halo. The value of the 
gravitational constant G is 1 in the standard units. 



In table 1 we summarized the model parameters and initial conditions of our iV-body 
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simulations. 

In our simulations, equations of motions of all particles in a dark halo and the satellite, 
i.e., N + 1 particles, are integrated self-consistently. In other words, the dynamical friction 
effect from halo particles to the satellite is included naturally. 

The number of halo particles N used in the simulations shown in this paper is 32768. 
We varied N from 8192 to 65536, and found any noticeable difference in the orbit of the 
satellite. 

We used GRAPE6 to calculate the acceleration. We adopted simple 0(N 2 ) direct 
summation, to avoid any possible numerical artifact caused by the approximations made in 
force calculation. BLLS used the treecode and JB used a composite grid-based code. We 
do not think the numerical method caused the difference, but we want to be absolutely sure 
that our iV-body simulation is as accurate as possible. 

We integrated the orbits of the satellites and halo particles using the standard leapfrog 
scheme with a shared stepsize, At = 0.03125. The error in total energy within 10~ 3 , which 
is small enough to observe the orbit evolution of satellites (Hashimoto et al, 2002). 



2.2. Semi-analytic Integration 

We performed semi-analytic calculations to follow the evolution of satellite orbits. 

In these calculations, the model of the satellite is the same as in the iV-body simulations, 
i.e., a single particle with mass M s and the softening length e s . 

Instead of being represented by N particles, the potential of the galactic halo is evaluated 
by using the gravitational potential of King 9 model with the same mass and scales as those 
adopted in iV-body simulations. 

In this integration, the force to the satellite due to the dynamical friction from the halo 
is evaluated by using an analytical formula, too. 

For the dynamical friction formula, we follow JB (and also Murai and Fujimoto) to use 
the standard "Chandrasekhar's dynamical friction formula" . It is expressed as 

^ = -l6n^m(M s + m) In A ^^ v, (3) 
at |v| d 

where m and M s are the masses of a component of host galaxy and its satellite galaxy 
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(Chandrasekhar 1943; Binney and Tremaine 1987). Here In A is the Coulomb logarithm 

lnA = ln( J R halo /e s \/ s 2 ), (4) 

where i?haio is the scale length of the galactic halo. This formula has been adopted by many 
semi-analytic studies of the orbital evolution of satellite galaxies (e.g., Murai and Fujimoto 
1980; Helmi and White, 1999; Johnston et al., 1995). It is also used in cosmological studies 
of galaxy formation in order to estimate the merging time scale of satellite galaxies (e.g., 
Kauffmann, et al, 1994). 

3. Result 

Figure 1 shows the orbital evolution of a model satellite galaxy. The ordinate and 
abscissa are the distance of the satellite from the center of the galaxy and time in the N- 
body units. The solid and dashed curves correspond to the result of A^-body simulation and 
that of semi-analytic model with standard dynamical friction formula (3). 

In Figure 1 two curves are in good agreement only for a first few dynamical times. 
After a few orbits, two curves deviate from each other. Figure 1 shows that the orbital 
decay calculated with formula (3) is faster than that obtained by A^-body simulation. If one 
measure the orbital eccentricity, it is clear that A^-body result shows only a small change in 
the eccentricity, while semi-analytic result shows significant circularization. 

Thus, even though we used completely different initial models and numerical method, 
we confirmed previous results by BLLS and JB that A^-body simulation shows little circu- 
larization while semi-analytical calculation with standard dynamical friction formula shows 
strong circularization. In the next section, we discuss the possible causes of this discrepancy. 

4. Possible causes of discrepancy 

Since we have obtained quite different results with A^-body and semi-analytic models, 
at least one of them must be wrong. Since A^-body calculation can suffer many numerical 
problems due to limited resolution and particle noise, one might think A r -body result is 
probably wrong. However, additional tests with different number of particles and different 
sizes of timestep showed very good agreement (Hashimoto et al, 2002). Therefore it seems 
our A^-body result is sound. In addition, as we stated in the previous section, our A^-body 
result is in good agreement with BLLS and JB. Though it is not impossible, it is certainly 
unlikely that all of these three works are wrong. 
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Fig. 1. — Time evolution of radius of satellite position from the galaxy center. Solid: result 
of iV-body simulation. Dashed: semi-analytical integration using constant A. 
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So let us now consider the possibility that the standard dynamical friction formula is 
wrong. 

The standard dynamical friction formula is obtained under the assumption that the 
massive object moves straight in a uniform and isotropic distribution of field particles. Field 
particles are also assumed to be moving straight, and any interaction between field particles 
is ignored. Clearly, the satellite does not move straight, but circle around the center of the 
parent galaxy. The distribution of field stars within the parent galaxy is far from uniform, 
and field stars also circle around in the parent galaxy. Thus, it is not really surprising that 
the naive use of the dynamical friction formula gives rather bad result. 

One obvious way to improve the accuracy of the dynamical friction formula is to calculate 
the linear response of the global distribution function of the parent galaxy to the presence 
and the orbit of the satellite (Weinberg, 1995). This approach would certainly give accurate 
and reliable result which agrees well with iV-body result (Hernquist and Weinberg 1989). 
However, since the global response depends on the distribution function itself, the result 
cannot be expressed in a compact and form. So here we consider the possibility to improve 
the standard formula. 

As we noted above, there are at least two problems with the standard formula. First, 
it assumes that both the satellite and field stars move straight. Second, it assumes that the 
density of the field star is the same everywhere. 

The first assumption is clearly wrong, but its effect is difficult to estimate. Let us 
consider the effect of the second assumption, which is much easier to evaluate. In previous 
works, the outer cutoff radius of the Coulomb logarithm is taken to be the scale length of 
the halo, while the representative density of the field stars is taken to be the local density 
around the satellite. This would clearly cause an overestimate of the Coulomb integral, for 
the case of the singular isothermal sphere (or the King model we used), since the stellar 
density drops off as fast as 1/r 2 . This means the logarithmic divergence of the Coulomb 
integral does not actually occur if we takes into account the effect of the density gradient. 

To correctly take into account the effect of the density gradient is a tricky problem. 
We cannot really use the straight line approximation for encounters with impact parameter 
comparable to or larger than R s , the distance to the center of the galaxy. On the other hand, 
it might not be too bad an assumption just to ignore any encounter with impact parameters 
comparable to or larger than R s . The density drops off rapidly and realistic effect is unlikely 
to enhance the effect of the encounter (except for the small fraction of the orbits in resonance 
with the orbit of satellite). 

Thus, it might be more sensible to use R s as the outer cutoff radius for the Coulomb 
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logarithm, that to use the traditional Rhaio- In fact, this use of R s is first proposed by a 
pioneering work by Tremaine (1976) on the effect of the dynamical friction to the orbit of 
LMC-SMC pair. 

To use e s as the inner cutoff is okay as an order-of-magnitude estimate, but can be 
improved by actually integrating the effect of all encounters with small impact parameters 
for Plummer potential, following the treatment by White (1976). For Plummer model, 
the integration can be performed analytically and the result is that inner cutoff radius is 
r in = 1.6e s . 

Figure 2 shows the result of the iV-body simulation compared to that of our improved 
(both outer and inner cutoffs) semi- analytic treatment. Figure 2 is the same as Figure 1 but 
for the choice of the Coulomb logarithm 



During our semi-analytic integration using (5), when the R s becomes smaller than 1.4e s , we 
simply put the dynamical friction term to be zero, since it is clearly unphysical to apply 
dynamical "acceleration" . 

Figure 2 shows that the results of the iV-body simulation and semi-analytic treatment 
agree quite well. 

The coefficient appeared in the denominator of the right-hand side of equation (5) is 
chosen in order to fit the curve of our semi-analytic model to that of ./V-body simulations. 
The Coulomb logarithm which we used to fit the result of direct iV-body simulation is in 
good agreement with the that proposed by White (1976). Strictly speaking, there is a slight 
difference between the value of the coefficient in the denominator we used in Figure 2 and that 
estimated by the method proposed in White (1976) : the former is 1.4, while the latter is 1.6. 
It is not serious. The latter value is estimated using straight-line approximation so that the 
value should be overestimated a little. (In other words, the straight-line approximation tends 
to underestimate the effect of dynamical friction.) We will discuss this problem elsewhere 
(Hashimoto, et al, 2002). 

The agreement between the iV-body result and semi-analytic treatment is quite remark- 
able. Figure 2 shows that the discrepancy shown in Figure 1 is caused by an inadequate 
estimate of A. Other possible reasons, such as the effect of the global response of the distribu- 
tion function, might still exists, but they are clearly not the prime reason of the discrepancy 
between iV-body and semi-analytic works which we discussed in the introduction and section 
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The improved agreement with the iV-body result is explained as follows. With b. 
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Fig. 2. — Same as Figure 1, but for the variable A. 
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R cu t, the semi-analytical treatment causes strong circularization and faster orbital evolution. 
This implies that the the semi-analytical treatment overestimated the dynamical friction 
around the perigalacticon. Around the apogalacticon, the error might exist, but relatively 
small compared to that at the perigalacticon. The use of variable b max reduces the value 
of In A both at perigalacticon and apogalacticon, but by a much larger factor at the peri- 
galacticon simply because R s is smaller. Thus, effectively we reduced the dynamical friction 
around the perigalacticon, which resulted in the improvement in the agreement with the 
iV-body result. 

In hindsight, it looks too obvious that the traditional use of the dynamical friction 
formula was inappropriate. Theoretically, it is clearly not justifiable to assume that the 
stellar density is the same up to the outer cutoff radius of the halo. From comparison 
between the A-body result and those of semi-analytic treatment, it also is clear that previous 
semi-analytic treatment overestimates deceleration due to the dynamical friction around the 
perigalacticon. 

To summarize our result, the orbital decay of satellites is slower than ever estimated, the 
eccentricity of orbit of revolution of a satellite around the host galaxy is almost constant. 
The reason why previous estimates are wrong is that previous studies overestimated the 
effect of dynamical friction at the perigalacticon. 



5. Summary and Applications 

We performed A-body simulations of satellite orbits. We found that the circularization 
of the orbit due to the dynamical friction is much slower than commonly believed. This 
discrepancy was also reported by BLLS, and we can see the same tendency from the numerical 
result reported by JB. 

Previous studies of satellite orbits used the outer cutoff radius of the dark halo as b max . 
We found that the effective b max should be of the order of R s , the distance of the satellite 
from the center of the galaxy, which varies as the satellite orbits around the galaxy. Our 
formula results in a greatly improved agreement with the A-body result. 



5.1. the Large Magellanic Cloud 

The Large Magellanic Cloud is the most famous satellite of Milky Way. Its orbit has 
been investigated from both observation and numerical simulations (e.g., Toomre, 1970; 
Tremaine, 1976; Lin and Lynden-Bell, 1977; Murai and Fujimoto, 1980). The importance 
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Fig. 3. — Radial Evolution of LMC. From -10 G years to 10 G years. 
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of the effect of dynamical friction from the galactic halo on the orbit evolution LMC is first 
emphasized by Tremaine (1976). 

By using numerical simulation, Murai and Fujimoto (1980) (hereafter MF) determined 
the orbital elements and the present phase of the LMC. They performed a number of back- 
ward numerical integrations of the orbits of the LMC and SMC from various initial con- 
ditions, and integrated orbits of test particles in the LMC and SMC for each condition. 
Comparing the result of distribution of test particles and the observed Magellanic stream, 
they chose the initial condition which gives the best fit. 

In their numerical integration, they assumed a halo expressed by a singular isothermal 
sphere, which is a simple flat-rotation halo. In their paper, it is not clear either what 
assumption or what exact value is adopted for In A, since there is no discussion on how they 
determined In A though it appeared in their equation (13). 

In order to see the effect of changing In A, we integrated the orbit of LMC both forward 
and backward in time, using both the constant A and variable A (b max = R s ). In this study, 
we express LMC as a single Plummer-softened particle with mass 2 x 1O 1O M and softening 
length 5 kpc. The rotation velocity of the halo is 250 km/s, same as what is used by MF. 
We simulated the orbit of the LMC only, since our purpose here is to demonstrate the effect 
of A and not the accurate determination of the orbits of the Clouds. 

The solid curve in Figure 3 corresponds to the orbit obtained when the dynamical 
friction is calculated using equation (5). The dashed curve Figure 3 correspond to the orbit 
obtained using the formula (3) and (4). Note that the backward part of this dashed curve is 
in very good agreement with the result of MF. This agreement strongly suggests that what 
MF used is indeed a constant A. 

Figure 3 shows that real evolution of the orbit of LMC (with variable A) is significantly 
smaller than what is obtained by MF. 10 Gyrs ago, the "true" apogalacticon was only 160 
kpc, while the solution by MF was 180 kpc. 

A more remarkable difference is in the future of the LMC. With the constant A. The 
LMC will fall to the galactic center in only 7 Gyrs with constant A, while our result suggests 
that it will take more than 14 G years for the LMC to fall to the galactic center. 

5.2. Statistical Evolution of Faint Galaxies 

In semi-anaritic studies of galaxy formation, it has been assumed that the orbits of 
satellite galaxies evolve through dynamical friction following Chandrasekhar's formula with 



-14- 



constant A. In this section, we discuss how our result might change our understanding of 
the statistical evolution of the satellite galaxies. 

Our study shows that the time evolution of the eccentricity of satellites is rather small. 
Thus, we may assume that the distribution of eccentricities of satellite galaxies at present 
directly reflects that at the formation epoch of the Galaxy. Therefore the distribution of 
eccentricities of satellites galaxies can be an important clue to the formation of the Galaxy. 

The lifetime of the satellite is estimated using the dynamical friction timescale with In A 
taken to be M H /M S (Lacey and Cole 1993; Kauffmann et al. 1994). This would cause a quite 
serious overestimate in the dynamical friction timescale, since the factor one should use is 
the ratio between the size of the halo and the size of the satellite. If we assume M oc a 4 , we 
have R oc M 1 / 2 . Thus, there is at least a factor of two difference in the value of In A. Since 
there are too many other uncertainties in the semi-analytic modeling of the galaxy number 
evolution, how serious this difference is not clear. However, it certainly affects the estimate of 
presently observed satellites rather strongly. A more detailed study on this aspect is clearly 
necessary. 
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thank Rainer Spurzem, T, Tsuchiya, and Andrea Just for helpful discussions. This work is 
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Table 1. Model Parameters of iV-body simulations 



Galactic Halo a M gal 


^halo c M sat 


esat 


Initial Position 


Initial Velocity 


King model (* = 9) 1.0 


0.0315 0.01 


0.1 


(1.5,0) 


(0,0.326) 



a Total mass of the host galaxy 

b Softening length of a halo particle 

c Satellite mass 

d Softening length of the satellite 
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Table 2. Model Parameters of LMC 



run Galactic Halo 


M ga i a V c M sat 


e aa t Initial Position 


Initial Velocity 


LMC Singular Isothermal Sphere 


6 7.3 x 10 n M o 250 km/s 2 x 1O 1O M 


5 kpc (50kpc ,0kpc) 


(50km/s,340km/s) 



a Constant circular speed of this galactic model 
b Halo mass within 50kpc from the galactic center 



